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Abstract 

In this paper, a method for the digital simulation of wind velocity fields by 
Fractional Spectral Moment function is proposed. It is shown that by constructing 
a digital filter whose coefficients are the fractional spectral moments, it is possi- 
ble to simulate samples of the target process as superposition of Riesz fractional 
derivatives of a Gaussian white noise processes. The key of this simulation tech- 
nique is the generalized Taylor expansion proposed by the authors. The method 
is extended to multivariate processes and practical issues on the implementation 
of the method are reported. 

1 Introduction 

Digital simulation of wind velocity field is needed in the design of wind exposed structures. 
The literature on this topic is vast and we refer to the review of Kareem, 2006 which provides 
a synthetic overview of different possible approaches. In this paper we confine our attention 
on the simulation of wind velocity fields which are both Gaussian and stationary. Under 
these assumptions, two strategies for the digital generation of wind velocity samples became 
standard: i) by superposition of trigonometric functions (Borgman, 1969; Shinozuka, 1971); 
ii) by digital filtering technique, which consists in calibrating the output of linear differential 
equation, called filter, excited by a white noise process. The latter are commonly referred 
to as auto-regressive (AR) algorithms, moving average (MA) and auto-regressive moving- 
average (ARMA) algorithms (see Deodatis and Shinozuka, 1988; Deodatis, 1995; Kozin, 1988; 
Naganuma et al, 1987; Saramas et al, 1985; Spanos and Mignolet 1986). 

This paper show a novel method for the representation of wind velocity fields by digital 
filtering. Yet, in contrast to the classical use of filters, that are linear differential equations, 
we propose to use fractional differential equation, in which fractional derivatives appears. 
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This is an absolute novelty in wind engineering and it will be shown that this method is 
easily applicable to multivariate fields. Our results rely on fractional calculus that is slowly 
becoming of interest in many engineering fields. For brevity's sake many concept on fractional 
operators will be given without further theory and readers are referred to the textbooks of 
Podlubny, 1998 and Samko et al. 1993 for introductive and authoritative treatments of such 
a topic. 

The method here proposed is based on the recently published papers of Cottone and 
Di Paola, 2010a and Cottone et al, 2010b, which posed the theoretical basis for a novel 
representation of stationary Gaussian processes as output of a fractional stochastic differential 
filter. 

In this paper, we extend our previous results both investigating on the digital simulation 
of samples with assigned statistics, useful in wind engineering, and applying the method to 
multivariate wind velocity field V (t) with assigned Power Spectral Density (PSD) matrix 
Sv (w). The paper is organized as follow. In section 2, fundamentals on the representation of 
multivariate wind velocity fields is given. Due to the Gaussian assumption, only second order 
statistics are needed to characterize the wind velocity field. This will be done by specification 
of the PSD matrix Sv (w), following the results given in Solari and Piccardo, 2001. Then, 
in sections 3-5, the representation method based on fractional spectral moments both for 
mono-variate and multivariate processes is given along with applications and implementation 
remarks, which aim to show the straightforward application of the theoretical concepts. 



2 Probabilistic Description of Multivariate Wind Ve- 
locity Fields 

In this section the wind velocity field characterization is recalled. Readers are referred to 
the paper of Solari and Piccardo, 2001 and references therein reported, for a comprehensive 
treatment of the topic. With the sake of introducing the notation, we consider that in the 
coordinate reference system x,y,z, the wind velocity process in a point P {x,y,z) can be 
expressed as 

(1) Vp{x,y,z-t) = V{z)+V{x,y,z-t) 

in which the mean value V (z), thought as function of the sole elevation, is added to the fluctua- 
tion around the mean component V {x, y, z; t) that is assumed to be a Gaussian stationary pro- 
cess with zero mean. The mean value follows the logarithmic profile V (z) = A;~^u^, In (z/zq), 
where k = 0.4 is the Karman's constant, is the shear velocity and zq is the roughness 
length. 

Given N points located in the space Pj {xj,yj,Zj) with j = 1,2, ...,N the wind velocity 
field can be represented as N-variate process by collecting the processes Vj {xj,yj, zj; t) in the 
vector V (t) = [FiV2...VAr] . Due to the Gaussian assumption, V (t) is characterized by the 
power spectral density (PSD) matrix Sv (w) that reads 



(2) 



Sv (t- 



SviVx (^) Sv^V2 {^) 

SviV2 (^) SV2V2 (^) 

SviVn {^) Sv2Vn {'^) 



-S'ViVat (^) 
SvnVn (^) 



In this paper, we will assume that only the co-spectrum, i.e. the real part of Sv {^), is needed 
to characterize the wind velocity field, neglecting the quad-spectrum, see Simiu &: Scalan, 
1996. The terms of the PSD matrix are calculated following Solari &: Piccardo, 2001, in 
which the hypotheses of flat homogeneous terrain and near neutral atmospheric conditions 
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are considered. Following the cited paper, explicit expression of the diagonal components of 
the PSD matrix reads 

c (,, dv alLv {z) /V {z) 



[l + 1.5dv ioLv {z) / {2irV {z))] 
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where Ly is the integral length scale of turbulence that defines the position of the turbulence 
spectral content, dy = 6.868, ay = I3u1 is the variance of the velocity longitudinal component 
and /? is a non-dimensional coefficient defined as turbulence intensity factor. 
The cross-power spectral components of Sv are expressed as 



(4) Sv^Vs {^) = \/SvrVr {^)SvsVs {^) exp {-frs (w)) 

where 



\ Cy iVr - Vsf + C| {Zr - Zg) 

(5) /,,(u;) = — V 



2^ {V{Zr) + V{Zs)) 

being Cy and Cz two coefficient to be experimentally determined. Parameters can be estimated 
accordingly to Solari & Piccardo, 2001. The method proposed in the next sections can also 
be applied to other form of the PSD matrix, but in this paper attention will be given to PSD 
constructed on eqs.(l3])-([5]). 

As the spectral decomposition of the PSD matrix will be needed in the next section, we 
will briefly introduce some notation here for reference's sake and some concept on the use 
of such a decomposition for the representation of wind velocity fields. Let us indicate with 
Ij {ijj) the eigenvalue associated with the eigenvector -i/'j (w) of the matrix Sv (w), * (w) is 
the eigenmatrix, in which the jth column is the eigenvector ipj [oj) and L {uj) is the diagonal 
matrix whose jth element is Ij {uj). Then, the relations 

(6) {uj) Sv (w) * (a;) = L (w) 
and 

(7) (w) * (cj) = I 

hold true, being I the identity matrix. The eigenmatrix decomposes the spectral density 
matrix such that 

(8) Sv (a;) = * (a;) L (a;) * (o;)'^ 

The spectral decomposition has been used by different authors, see Di Paola, 1998 and Di 
Paola and Gullo, 2001, to express the wind field N-variate processes V (t) in the form 

N oo , 

(9) V (t) = ^ / (u) J I, (u) e^-'dB, (u) 

where Bj (u) is a zero mean normal complex process having orthogonal increments, i.e. 
E[dBj{uj)] = 0, dBj{uj) = dBj (to) and E [dBj (uJr) dB"^ (uJs)] = ^i^ri^s^jkd^r-, having indi- 
cated by the over bar the complex conjugate and by 5jk the Kronecker delta (i.e. 5jk = if 
i i^k and = 1 if j = fc). By eq.(l9]) the vector process V (t) is decomposed into the sum of 
coherent and independent elementary vectors 



(10) Yj{t)= (w) Jlj (w) e*"*dBj (w) 

J —CXD 

such that 

N 

(11) Y{t) = Y,Y,{t) 
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Moreover, retaining only the first M <^ N most significant eigenvalues the latter is approxi- 
mated in the form 

M 

(12) v(t)-J]y,(t) 

As pointed out by Di Paola, 1998, although the spectral decomposition implies the calculation 
of the frequency dependent eigen-properties of the power spectral matrix Sv(w), the com- 
putational pay off is that only M <^ N independent component Yj can be considered. This 
concept can be combined also to the method here proposed as shown in the following. 

3 Fractional Spectral Moments Decomposition 

In this section we start considering a Gaussian stationary mono-variate stochastic process 
with known statistics, while in the next section the proposed method is extended to Gaussian 
stationary multi-variate processes. Let us assume we want to simulate a sample of the wind 
velocity in a point, with assigned PSD Sy (w) whose form is given in the previous section. The 
method here proposed follows two steps. Firstly, we will express an assigned PSD function in 
terms of fractional spectral moments. Then, we will construct a digital filter whose output 
has the assigned PSD. It will be shown that fractional spectral moments of the linear system 
transfer function are the filter coefficients and that the filter has the form of a fractional 
differential equation. The hinge of the method is the generalized Taylor expansion introduced 
by the authors and applied in different contexts such as in probability, in Cottone and Di 
Paola, 2008, in stochastic dynamics by path integral solution in Cottone et al. 2009, in the 
solution of stochastic differential equations in Cottone et al. 2008. The application of the 
generalized Taylor expansion to representation of power spectral densities and correlation 
function has been introduced in Cottone and Di Paola, 2010, and Cottone et al. 2010 in the 
mono-variate case and readers are referred to those papers for more insight on the topic and 
relevant demonstrations. 

3.1 Fractional moments for stochastic processes 

Let us consider a stationary Gaussian stochastic process V{t) with target Sv{i^)- Recall the 
definition of the Spectral Moments (SMs) of the process, see Vanmarcke, 1972, defined as the 
integral 

poo 

(13) \{r = / oj'JGv (w) doj j = 0, 1, 2, ... 

Jo 

where Gy (w) is the one-sided PSD calculated as Gv{u}) = 25y(ci;)?7(a;), being U{ijj) the unit 
step function. Physical meaning of SMs quantities has been given in Di Paola (1985). In 
Cottone and Di Paola, 2010a, the authors introduced a generalized class of SM, called the 
Fractional Spectral Moments (FSMs) of the one-sided PSD, that is 

POO 

(14) Ay (7) = / uj^Gv (oj) duj 7 E C 

Jo 

Although the use of the adjective "fractional" to indicate moments that are of complex order 
7 might be confusing, we keep it for similarity with the consolidated terminology of the 
Fractional Calculus which represent the calculus for derivation of real or complex order. The 
FSM function Ay (7) is a complex function that keeps all the information to restores both the 
power spectral density Sy (w) and the correlation function Ry (r). Skipping the details on 
the demonstration, see the paper above cited, we report here just the final formula needed in 
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the following. To this aim, let us introduce the Riesz fractional integral and derivative as 
21^(7) J-00 It -si 



(15a) (r/)(i) = __/ iie7> 0,7/1,3,5,... 



(15b) (Vy) it) = r li^-^l-JM^s Re^ > 0, 7 / 1, 3, 5, ... 

where (7) = F (7) cos being F (7) the Euler gamma function, ^ = p + ir], (7, r/ G M), 

and i = \/— T. 

It has been shown that the PSD and the correlation function are expressed in terms of 
FSMs as 

1 rp+ioo 

(16) Ry(r) = — V (7) ky (-7) |r|-T d7 

^VTi J p—ioo 

and 

1 rp+ioo 

(17) Sv{u) = — Ay(-7)|a;r-'d7 

J p—ioo 

Both integrals are valid for p belonging to an interval that depends on the convergence 
conditions of eq. ()14p . which simply becomes the interval < p < 1 for absolute integrable 
functions. These equations are a form of the generalized Taylor expansion for functions which 
are both symmetric and Fourier pairs. Not-symmetric functions are treated in the paper 
Cottone and Di Paola, 2009a. Eqs. ()16p and ()17p are integrals along an axis that is parallel 
to the imaginary axis, called Bromwich's path, and can be simply approximated. Indeed, it 
suffices to truncate the integration in the interval [—r]s,i]s] and to evaluate the integrand in 
2m + 1 points located at steps Arjs = ijs/m. Then, posing 7^, = p + ikAij, with k = — m, m, 
the approximated forms of eqs. (fT6]l and (fT7|l are 

A 

(18) ^i^k)Av{-7k)\rr'' 

k=—m 



(19) Sv{oj) = ^ Yl Av{-7k)m^''-' 



and 

At? 
47r 

k=—m 

The following simple example clarifies the application of the previous formula. 
3.1.1 Example 

Let us consider a simpler form of eq.([3]) 
(20) Sv (w) 



[l + b\uj 



15/3 



with a > and 6 > 0. By using Wolfram's Mathematica software it is easy to calculate FSMs 
by applying the definition reported in eq. p^ . thus obtaining 

. . ^ 2a6-(^+^)F(2/3-7)F(l + 7) 

only if —1 < p < 2/3, having care that in the whole paper p indicates the real part of the 
complex variable 7. As in eqs. ()16p and ()17p we use Ay (—7), the fundamental strip in which 
such equations are valid is — 2/3 < p < 1. Inside this interval, we can choose any value of p in 
order to reconstruct the functions Sy (aj)and Ry (t), because the integrands are holomorphic. 
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Figure 1: Comparison between the exact (continuous line) and the approximated (dot- 
ted hne) correlation function, left panel in log-log diagram, and between the exact 
(continuous line) and the approximated (dotted line) PSD, right panel 



Outside this interval, the residue theorem must be applied. This choice can be profitable and 
it has been used in Cottone and Di Paola, 2010a to regularize the approximated formula in 
eqs. pH]) and (HHD in zero. 

Figures 1 shows the application of eq. ()18p and ()19p with a = 374.8 and b = 4.51 which 
corresponds to having chosen the parameters zq = 0.7, = 2m/s, /3 = 4.96,(7^ = 10, 
z = 5m and Ly (z) = 300 (z/200)°-^^+°-°^'''(^°Hn eq.Q. Sixty FSMs of the type Ay {-Jk) = 
1/2+ik (0.1) with k = —30, 30 have been calculated and stored from eq. ()2ip and introduced 
in eqs.([18]) and (1191) . 

This simple example shows that a finite numbed of FSMs can be used to represent the 
statistics of the Gaussian stationary process V{t). On this solid ground, we can simulate the 
process V{t) by a fractional linear filter whose coefficients are FSMs. 

4 Simulation of stationary Gaussian stochastic pro- 
cesses by FSMs 

Objective of this section is to represent a Gaussian stationary process V{t) with assigned PSD 
(S'y((^) as the output of a fractional differential equation. As shown in Cottone et al. 2010, to 
this aim, it suffices to consider an ideal linear system C{V{t)) = W{t) where £(•) is a linear 
operator and W{t) is a Gaussian white noise process with zero mean and correlation function 
E[W(t)W{t + r)] = q5{t), and power spectral density Sw = Q'/(2'/r), being q the intensity 
parameter. From linear system theory, it is known that the output V{t) can be characterized 
both by the impulse response function h{t) through the Duhamel integral 

(22) V{t)= [ h{t-T)W (r) dr 



or by its Fourier transform H[uj), namely the transfer function, through the input-output 
relation that reads 

(23) Sv{uj) = \H{uj)\^Sw{u^) = ^\H{u)\^ 

Assuming Arg [H [uj]) = 0, we get a non causal differential equation characterized by the 
transfer function 

(24) H{io) = \H{io)\ = ^'^-^Sv{u) 
We calculate now the fractional spectral moments of H{uj) as 

(25) IiH{-l)'^= f H{uj)Auj, Re7>0 
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that are called H-Fractional Spectral Moments (H-FSMs). As before, it can be shown that 
such a function can be used to represent both the impulse response in the time domain as 

(26) h {t) = -— ^ / u (7) (-7) t-^dj, t > 
and the transfer function in the form 

1 j-p+ioo 

(27) H{io) = — nH(-7)l^r-'d7 

^■^^ J p—ioo 

Following the same notation of the previous section, approximated form of the last integrals 



are 

(28) ^ (*) - TT^ E ^ (^^) (-^fc) 



TT / M |7fe-l 



(29) ^(^)= 4^ E ^Hi-lk) 

k=—m 



As shown in Cottone et al. 2010, the input-output relation for linear system and Fourier 
transform properties of Riesz fractional operators lead to the relevant representation of the 
Gaussian stationary process 



1 pp+ioo 

(30) v{t) = — Uh (-7) {I'~^W) (t) d7 

"^^^ J p~ioo 



f p+ioo 
I p—ioo 

whose PSD is Sv{oj)- The approximated form reads 



(31) v{t) = ^Y.^^ (-^'^) (t) 

k=~m 

and readers must bear in mind that this approximation carries out a truncation and discretiza- 
tion error, that can be made arbitrarily small. Eq. ()30p . or its discretized counterpart given 
in eq. pT]) . show that the process V{t) may be obtained as the superposition of the fractional 
integrals of a Gaussian white noise process. 

4.1 Implementation 

We continue the previous example by simulating a digital temporal signal by means of eq. (|31|) . 
FSMs have been already calculated and, for the functional form in eq. ()20p they are given by 
eq. (|21|) . It is reasonable to keep the parameters p = 1/2, Ar; = 0.1 and m = 30 by which, in 
the previous example, we achieved to represent both the PSD and the correlation function, 
to a good accuracy. Firstly, we must calculate the H-FSMs, and this task can be strongly 
simplified by using of the Wolfram's Mathematica, which returns 

2^/2^6"(i+^)r(-l/6-7)r(l + 7) 

(32) lin (7) = 

Then, the practical implementation of eq. ()3ip requires a further step to calculate the Riesz 
fractional integral of the Gaussian white noise W {t). To this aim, we firstly exploit the 
equivalence between the Riesz's fractional integral and the Grnwald - Letnikov's discrete 
fractional operator (see Samko et al. 1993) that in approximated form is expressed as 
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where 

(34) {AlW) (t) = f; (-l)^' (^l^WitTkr) 

being r a discrete time step. To the hmit r — )• 0, eq. ()33p becomes an identity. These relations 
are true on the whole real support and they can be also rewritten in the case of a finite interval. 
Let us consider a time window [0, T] partitioned into n time steps of amplitude r. At each 
node 0, 1, n, the noise W {0) ,W Qi) , W (jr) , W (nr) is evaluated as a realization of a 
Gaussian random variable with zero mean and variance \/r, having chosen g = 1. The Riesz's 
fractional integral of such a noise is calculated by the Grnwald - Letnikov's approach as 

W {jr -kT)+Yl (-1)' [k^ )^^^^ + 
(35) 

The coefficients in the latter can be efficiently calculated both iteratively and by using the Fast 
Fourier transform as reported in Podlubny, p. 209. Such coefficients, calculated for Re7 > — 1, 
decrease with inverse power-law behavior and, for many functions, they can be neglected after 
a finite number of terms, say p. In the case here considered, with p = 1/2, p = 400 terms 
suffice to achieve the searched accuracy. In other words, the fractional integral of the Gaussian 
white noise at the time step jr is calculated by eq. ([35|) considering the influence of p noise 
realizations which precede and follow the instant j. 




5 10 15 20 25 30 

time lag r [sec] 

Figure 2: Comparison between the exact (continuous line) and the simulated (dotted 
line) correlation function. 



(rW) (jr) 



2cos (77r/2) 



Figure 2 shows the comparison between the target correlation function and the correlation 
function of the time series V {jr) with j = 0, ...,n. The simulation has been implemented in 
Octave, see Eaton, 2002, with the parameters: 3 • 10^ time steps and r = 5 • 10"^. 



5 Extension to multivariate stochastic processes 

Let us consider a vector of Gaussian white noise processes, W(t), with zero mean and diagonal 
power spectral density matrix Sw('^)) whose elements along the diagonal are SwjWj ('^) = 
Qj/ (27r) , where qj is the strength of the jth noise component. We can model V (t) as the output 
of a multidimensional linear system characterized by the impulse response matrix H(t) that 
can be expressed by Duhamel's integral 

(36) Y{t)= f H(t-s)W(s)ds 

J — oo 
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Indicating by Sv (i^) the power spectral density matrix of the output, from the hnearity of 
the system it follows that 



(37) Sv (lo) = H (u) S^' (oj) S^' (oj) (oj) 

where H (cj) is the transform matrix and the bar means conjugation. Rewriting eq.(l6]) as 

(38) Sv (uj) = * (u) (^) Li/2 

and comparing the latter with eq. ()37p the characterization of the linear system can be obtained 
by selecting the transfer function as 

(39) U{uj) = ^{uj)L^/^{uj) 

with the choice of qj = 27r for every noise component Wj (t).The output of course will be 
Gaussian due to the linearity of the system. 

The next steps we want to take aims to represent V (t) in terms of fractional spectral 
moments of the transfer matrix. This will turn out to have relevant computational advantages. 
Firstly we will introduce the fractional moments of the transfer matrix H(a;); then we will 
represent H(w) as a sum of fractional spectral moments and finally we will give the expression 
of the process V (t) with assigned Sv (w) in terms of the fractional spectral moments of H(a;). 

Let us define the H-FSMs of H(tj) of complex order j = p + ir], labeled in the following 
by n(7), as 

(40) n(7) =^ /"°° |wrH(a;)dw 

J — oo 

with p chosen such that the integral converges. 

In perfect analogy to the eqs. (f28]) and ([29]) we can obtain the representation of the impulse 
response matrix H(t) in the time domain as 

-j^ rp+ioo 
I p—ioo 



(41) H (t) = — / 1/ (7) n (-7) t-Mj 



and the transfer function H(ti;) in the form 

(42) H(a;) = — / Il{-^)\u\^-U-f 

Assuming that all the components of the transfer matrix H (w) are absolute integrable, in 
the following we will consider the fundamental strip of eqs. (|41|) and (|42p such as < p < 1. 

A simple approximation of the latter relations can be given by truncating the integral 
along the imaginary axis r]. Indicating as rjs the truncation limit for the integrals, let us 
divide the interval [—i]s,ris] in 2m intervals of amplitude A77 = r]s/m, with m G N. Then, the 
integrals in eqs. (j^T]) and (P2|) can be approximated by the values at the nodes 'Jk = P + ikArj 
in the form 

(43) Hrs{t)^^ Yl z^(7fc)n.s(-7fc)t-^^ 



(2^)^ 



-m 



and 

(44) Hrs (^) = ^ E ^rs i-lk) H^'^' 

k=—m 

The latter equations hold true for each component r, s = 1, N of the matrices H(t) and 
H(^). 
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Now, having represented the transfer function both in exact and approximated form in 
terms of H-FSM we are ready to infer an analytic expression for the representation of the 
stationary Gaussian vector process V(t). 

First, the input-output relation of eq. ()32p is written in Fourier domain as \^ {uj,T) = 
H (cj) W (w, T), where T > is a truncation bound. Then, by introducing the representation 
of the transfer matrix in terms of H-FSMs, one obtains 

1 rp+ico 

(45) V{oo,T) = — n(-7) W(w,T)d7 

4vr? J p-ioo 

and by inverse Fourier transform, we finelly obtain the relation searched that reads 

1 r-p+ioo 

(46) V (t) = — / n (-7) (/1-^w) {t) 

In the previous derivation we have used the property of the inverse Fourier transform of 
the Riesz's fractional integrals 



lim F-^ \\uj\^~^ Wj {io,T);t] = {l^'^'WA (t) 



that is true for each component of W (t) . 

Eq. (|46p is a novel exact representation of the multivariate stationary Gaussian process 
with assigned PSD matrix, in which the term 

(47) {I^'^W) (t) = [ I^-^Wi I^-^W2 ... I^-^Wn f 

is a vector whose elements are the Riesz fractional integrals of complex order 1 — 7 of the 
independent white noise processes Wj, j = 1, 2, N. The fundamental strip in which Eq. ()46p 
converges, depends on the integrability of the single components of Eq. (j40p . Following our 
assumption that the components of the matrix H (a;) are absolute integrable, the integral in 
Eq. ()46p can be performed inside the strip < p < 1. As already pointed out, any values of p 
inside this fundamental strip can be chosen. 



5.1 Approximated decomposition by H-FSM 

Eq. ()46p is approximated in a quite simple but effective way, following the same trace of the 
approximations in eqs. ()43p and (|44p . To this aim, a value of p is chosen inside the interval 
< p < 1. Then, the integrator is rewritten as = idr] and the bounds of integration are 
truncated up to the values [— r/s,f/s]. In this way we are introducing a truncation error in the 
evaluation of V(t), and the integral becomes 

(48) r " ^"^^ (^'""W) it) dri 

The interval \—r]s,r]^ is partitioned in 2m intervals of amplitude A?? = r]s/ra, with m G N 
and the integral in eq. fj^H]) is evaluated at the nodes 7^ = p + i/cAry. This correspond to apply 
a rectangular numerical scheme to the evaluation of eq. (j35|) . that reads 

(49) V (i) - ^ n (-7.) (/^-^^ W) (t) 

fc=— m 

The accuracy of the latter formula can be increased by using more accurate numerical schemes 
to approximate the exact decomposition in eq. (|46p . This kind of problem is known in the field 
of the numerical treatment of the Laplace transform inversion, that has the same mathematical 
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structure of the problem here treated. Talbot's contours, conformal mapping, higher numerical 
integration schemes are of course applicable to our problem, but the simplicity and accuracy of 
eg . (1491) are satisfactory and no further numerical effort is needed to our scopes. It is important 
to emphasize that to practically use eq. ()49p . (2m+l) x A^^numbers, i.e. the H-FSM at different 
values of fc = — m, ...,m, must be calculated. This number can be strongly reduced by the 
spectral decomposition of the PSD matrix recalled in section 2. 

5.2 H-FSM in reduced space 

The computational effort for the calculation of the H-FSM decomposition in eq. ()46p is mainly 
influenced by the calculation of the fractional spectral moments. A considerable reduction of 
computational effort can be achieved by considering the eigen-properties of the PSD matrix. 
Indeed, following the paper of Di Paola and GuUo, 2001, in many applications, like in wind 
engineering, only a reduced number of terms M ^ is relevant in simulating the N-variate 
process V(t), in the form of eg . (fT2]l . 

Applying this reasoning to the H-FSM decomposition, a reduced form of the integral 
representation by H-FSMs can be given in the form 

1 fp+ioo 

(50) V (i) = — / n (-7) (/1-Tw) (t) d-f 

where H (—7) is a (A^ x M) reduced matrix and ^/^^'''W^ {t) is a (M x 1) vector and fl (—7) 
is calculated on the reduced matrix H (w) as 

By such a reduction, the discrete approximated form of the latter reads 

(51) V (t) - ^ J] n (-7.) [I'-^-w) it) 

k=-~m 

and the computation is performed storing (2m + 1) x x M numbers, with M <^ N. 

6 Conclusions 

We have introduced a novel method for the representation and the consequent digital sim- 
ulation of stationary Gaussian processes and multivariate fields. The methods is developed 
introducing the Fractional Spectral Moments. If these features are calculated from the power 
spectral density function, then they can be used to represent both in exact and in approximate 
form the correlation and the spectral density, as shown in section 3. If the Fractional Spectral 
Moments are calculated by integrating the transfer matrix of a linear system excited by a 
Gaussian white noise, then they are the coefficients of the time series given in eq. ()3ip which 
restore the target process. This method can be seen as a valid alternative to classical ARMA 
models, especially in case of PSD function with pathological behavior, see Cottone and Di 
Paola, 2010a. 

Moreover the paper extends the use of the fractional spectral moments to represent N- 
variate Gaussian processes. Two steps are required to this aim. The first step consists in 
performing a spectral decomposition of the PSD matrix of the N-variate process, retaining 
only the most relevant Ad <ti N eigenvectors. Then, such quantities are used to calibrate 
the transfer matrix of a reduced multidimensional linear system whose output is the searched 
process and the input are M uncorrelated Gaussian white noise processes with unitary power 
spectral densities. 
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